#!/bin/bash
#从 kingdom 导出层位然后 line drawing,由于 horizon 需要新建很多层位,因此用断层表示层位,导出的都是断层格式(GeoQuest-CARD IMAGE 7),注意左对齐
function preset()
{
    #=================================
    it=0
    line=04
    fig_name=${line}_bend_fault
    fig_fmt=png
    dmig=../../sgy/grd/${line}_migkd_it${it}.grd
    #======bend fault all (used for connetc dashed line)=========
    echo 208722.44 6131 6329 >tmp.txt
    echo 211233.89 6238 6610 >>tmp.txt
    echo 216196.46 6375 6465 >>tmp.txt
    echo 219370.43 6672 6738 >>tmp.txt
    echo 223268.67 6853 7101 >>tmp.txt
    echo 226317.83 6993 7092 >>tmp.txt
    echo 230769.20 7216 7319 >>tmp.txt
    echo 232844.53 7443 7616 >>tmp.txt
    echo 234405.85 7588 7674 >>tmp.txt
    echo 236588.67 7612 7810 >>tmp.txt
    echo 241896.79 8462 8631 >>tmp.txt
    echo 245521.33 9089 9242 >>tmp.txt
    echo 248805.93 9448 9626 >>tmp.txt
    awk '{print$1, $3-$2}' tmp.txt >${line}_bend_all.txt
    #=============bend fault_oceanic_landward====================
    echo 208722.44 6131 6329 >tmp.txt
    echo 211233.89 6238 6610 >>tmp.txt
    echo 216196.46 6375 6465 >>tmp.txt
    echo 219370.43 6672 6738 >>tmp.txt
    echo 223268.67 6853 7101 >>tmp.txt
    echo 226317.83 6993 7092 >>tmp.txt
    echo 230769.20 7216 7319 >>tmp.txt
    awk '{print$1, $3-$2}' tmp.txt >${line}_oceanic_landward.txt
    #=============bend fault_oceanic_seaward====================
    #=============bend fault_underplate_landward====================
    echo 232844.53 7443 7616 >tmp.txt
    awk '{print$1, $3-$2}' tmp.txt >${line}_plate_landward.txt
    #=============bend fault_underplate_seadward====================
    echo 234405.85 7588 7674 >tmp.txt
    echo 236588.67 7612 7810 >>tmp.txt
    echo 241896.79 8462 8631 >>tmp.txt
    echo 245521.33 9089 9242 >>tmp.txt
    echo 248805.93 9448 9626 >>tmp.txt
    awk '{print$1, $3-$2}' tmp.txt >${line}_plate_seaward.txt




}

preset

gmt begin $fig_name $fig_fmt E300 
# ============================================================ (a) seismic profile===========================================================================================
origin_x_min=2630 #set min_cmp num 
origin_x_max=15500 #set max_cmp num 
new_x_min=`echo "scale=2; $origin_x_min * 6.25 / 1000 " | bc`
new_x_max=`echo "scale=2; $origin_x_max * 6.25 / 1000 " | bc`
origin_d_min=0
origin_d_max=20
xylimits_cmp=`echo $origin_x_min/$origin_x_max/$origin_d_min/$origin_d_max`
trench_cmp=6900
trench_km=`echo "scale=2; $trench_cmp * 6.25 / 1000 " | bc`
left_trench=`echo "scale=2; $new_x_min - $trench_km " | bc`
right_trench=`echo "scale=2; $new_x_max - $trench_km " | bc`
xylimits_km=`echo $left_trench/$right_trench/$origin_d_min/$origin_d_max` #6900*6.25/1000
height=`echo "scale=2; ($origin_d_max - $origin_d_min)/3" | bc`
vert_exa=1
width=`echo "scale=2; $height * (${origin_x_max} - ${origin_x_min}) * 6.25 *0.001  / ((${origin_d_max} - ${origin_d_min})* ${vert_exa})" | bc`
height_new=`echo "scale=2; ${height}/3" | bc` #used for narrow height but the width is same as profile length, this is used for the offset of bending fault
throw_min=0
throw_max=1000
y_shift=`echo "scale=2; $width+0.3 " | bc`
x_shift=`echo "scale=2; $height+0.8 " | bc`
for v in ${line}  #读出起始终止坐标
do
x1=$(awk '{if($1=='$origin_x_min') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
y1=$(awk '{if($1=='$origin_x_min') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
x2=$(awk '{if($1=='$origin_x_max') print $2}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
y2=$(awk '{if($1=='$origin_x_max') print $3}' ../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
done

gmt basemap -R$x1/$x2/${throw_min}/${throw_max} -JX${width}c/${height_new}c -Byf50a200+l"Fault throw (m)" -BWesn  \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p

#=============all bend====================
gmt plot ${line}_bend_all.txt -W0.03c,grey,- 
if test -f "${line}_oceanic_landward.txt"; then
gmt plot ${line}_oceanic_landward.txt -W0.1c,blue -Ss-0.1c
fi
if test -f "${line}_oceanic_seaward.txt"; then
gmt plot ${line}_oceanic_seaward.txt -W0.1c,blue -St-0.1c
fi
if test -f "${line}_plate_landward.txt"; then
gmt plot ${line}_plate_landward.txt -W0.1c,cyan -Ss-0.1c
fi
if test -f "${line}_plate_seaward.txt"; then
gmt plot ${line}_plate_seaward.txt -W0.1c,cyan -St-0.1c
fi
gmt basemap -R$xylimits_km -JX${width}c/-${height_new}c -Bxf1a5+l"Distance to trench (km)" -BS \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt basemap -R$xylimits_cmp -JX${width}c/-${height_new}c -Bxf100a1000+l"CMP number" -BN \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt text -F+cTC+f5p,1 -Dj0.1c/0.1c <<EOF
MC${line}
EOF
gmt text -F+cTL+f5p,1 -Dj0.1c/0.1c <<EOF
W
EOF
gmt text -F+cTR+f5p,1 -Dj0.1c/0.1c <<EOF
E
EOF


gmt end show
rm tmp*
rm gmt*